setwd("")
library(foreign)
library(dplyr)
source("functions.R")


################################## ################################## 
################################## ################################## 
############ Figures
################################## ################################## 
################################## ################################## 


################################## ##################################  
################################# Figure 1: Correlational Plots

## Note: To reduce outliers, small municipalities (<500 pop), as well as municipalities at the upper tail of immigrant settlement are removed (2.5% tail)

pdf("Figures/Figure1.pdf",width=8, height=5)
par(mfrow=c(2,5),mar=c(3,2,2,.8))
loessplot("./Datasets/Crossnational/Belgium.dta",0,"Belgium",.975,20)
loessplot("./Datasets/Crossnational/Denmark.dta",0,"Denmark",.975,35)
loessplot("./Datasets/Crossnational/Netherlands.dta",0,"Netherlands",.975,30)
loessplot("./Datasets/Crossnational/Norway.dta",0,"Norway",.975,20)
loessplot("./Datasets/Crossnational/Sweden.dta",0,"Sweden",.975,35)
loessplot("./Datasets/Crossnational/France.dta",1,"France",.975,15)
loessplot("./Datasets/Crossnational/Italy.dta",1,"Italy",.975,25)
loessplot("./Datasets/Crossnational/Spain.dta",1,"Spain",.975,20)
loessplot("./Datasets/Crossnational/Swiss_German.dta",1,"Switzerland",.975,60)
dev.off()


################################## ##################################  
################################# Figure 2

pdf("Figures/Figure2.pdf",width=8, height=4)

par(mfrow=c(1,2))
data<- read.dta("./Datasets/Other/Swiss_master.dta")
# Due to high variance in municipal size in Switz, small population trim to avoid distorting growth rates at edges of distribution
data <- subset(data, data$pop >= quantile(data$pop,.05,na.rm=T) & pop <= quantile(data$pop,.95,na.rm=T)) 
output <- matrix(NA,6,5)
data$agroup <- data$foreignbin4
dv <- data$net_welfare_head / (data$priceindex/100)

output[1,] <- c(1,tapply(dv[data$reformyear==-3],data$agroup[data$reformyear==-3],mean,na.rm=T))
output[2,] <- c(2,tapply(dv[data$reformyear==-2],data$agroup[data$reformyear==-2],mean,na.rm=T))
output[3,] <- c(3,tapply(dv[data$reformyear==-1],data$agroup[data$reformyear==-1],mean,na.rm=T))
output[4,] <- c(4,tapply(dv[data$reformyear==0],data$agroup[data$reformyear==0],mean,na.rm=T))
output[5,] <- c(5,tapply(dv[data$reformyear==1],data$agroup[data$reformyear==1],mean,na.rm=T))
output[6,] <- c(6,tapply(dv[data$reformyear==2],data$agroup[data$reformyear==2],mean,na.rm=T))
ulim <- max(output)
llim <- min(output[,2:6])

output2 <- output
scaleyear <- 1
output2[,2] <- (output2[,2] / output2[scaleyear,2]) * 100
output2[,3] <- (output2[,3] / output2[scaleyear,3]) * 100
output2[,4] <- (output2[,4] / output2[scaleyear,4]) * 100
output2[,5] <- (output2[,5] / output2[scaleyear,5]) * 100
output2
plot(output2[,1],output2[,2], type="n",lty=3,ylim=c(98,140),main="Switzerland",xlab="",xaxt="n",ylab="Average Expenditure (Indexed)")
grid(col="gray92")
rect(3,0,4,200,col="gray95",border="NA")
axis(1, at=output2[,1],tick=FALSE,hadj=.1,labels=c("t-3","t-2","t-1","t","t+1","t+2"),las=1,cex.axis=.88)
lines(output2[,1],output2[,2],lty=2)
box()
lines(output2[,1],output2[,5],lty=1)
legend("topleft", c("Lowest Qtile", "Highest Qtile"), col=c("black","black"),lty=c(2,1),cex=.7,bg="white")
box()


## Belgium
belgiumdata<- read.dta("./Datasets/Other/Belgium_master.dta")
output <- matrix(NA,6,5)
z <-1
dv <- belgiumdata$cpas_head
for (i in  2002:2007){
  output[z,] <- c(i,mean(dv[belgiumdata$year==i & belgiumdata$foreignbin4==1],na.rm=T),
                  mean(dv[belgiumdata$year==i & belgiumdata$foreignbin4==2],na.rm=T),
                  mean(dv[belgiumdata$year==i & belgiumdata$foreignbin4==3],na.rm=T),
                  mean(dv[belgiumdata$year==i & belgiumdata$foreignbin4==4],na.rm=T)
  )
  z <- z+1
}

output2 <- output
scaleyear <- 1
output2[,2] <- (output2[,2] / output2[scaleyear,2]) * 100
output2[,3] <- (output2[,3] / output2[scaleyear,3]) * 100
output2[,4] <- (output2[,4] / output2[scaleyear,4]) * 100
output2[,5] <- (output2[,5] / output2[scaleyear,5]) * 100
output2
plot(output2[,1],output2[,2], type="n",lty=3,ylim=c(100,111),main="Belgium",xlab="",xaxt="n",ylab="Average Expenditure (Indexed)")
grid(col="gray92")
rect(2004,0,2005,400,col="gray95",border="NA")
axis(1, at=output2[,1],tick=FALSE,hadj=.1,labels=c("t-3","t-2","t-1","t","t+1","t+2"),las=1,cex.axis=.88)
lines(output2[,1],output2[,5],lty=1,col="black")
lines(output2[,1],output2[,2],lty=2,col="black")
legend("topleft", c("Lowest Qtile", "Highest Qtile"), col=c("black","black"),lty=c(2,1),cex=.7,bg="white")
box()
dev.off()



################################## ################################## ################################## 
################################## Figure 3: Politics

pdf("Figures/Figure3.pdf",width=8, height=4)
par(mfrow=c(1,2))

swissdata <- read.dta("./Datasets/Other/Swiss_master.dta")
swissdata <- subset(swissdata,swissdata$reformyear==1)   
swissdata <- subset(swissdata,swissdata$interaction_constant >= quantile(swissdata$interaction_constant, c(.33),na.rm=T)) # Limit to top 2 terciles
swissdata <- subset(swissdata,swissdata$pop > 250) # remove hamlets
swissdata <- subset(swissdata, twoyear_growth > quantile(swissdata$twoyear_growth,.01,na.rm=T) & twoyear_growth < quantile(swissdata$twoyear_growth,.99,na.rm=T))  # Outlier trim
mdata <- na.omit(data.frame(swissdata$cantonid,swissdata$left,swissdata$twoyear_growth))
x <- mdata$swissdata.left
y <- mdata$swissdata.twoyear_growth 
y <- y* 100
general_loess(x,y,1.5,xlab="% Left-wing Vote Share",ylab="% Change in Social Spending",c(-5,25),"Switzerland",0,.975)
abline(h=0)

belgiumdata <- read.dta("./Datasets/Other/Belgium_master.dta")
belgiumdata <- subset(belgiumdata,belgiumdata$year==2005)
belgiumdata <- subset(belgiumdata, belgiumdata$main_interaction >= quantile(belgiumdata$main_interaction, c(.33),na.rm=T))  # Limit to top 2 terciles
mdata <- na.omit(data.frame(belgiumdata$leftshare_begin,belgiumdata$twoyear_growth))
x <- (mdata$belgiumdata.leftshare_begin*100)
y <- mdata$belgiumdata.twoyear_growth * 100
general_loess(x,y,1.5,xlab="% Left-wing Council Seats",ylab="% Change in Social Spending",c(-5,25),"Belgium",0,.975)
abline(h=0)
dev.off()

################################## ################################## ################################## 
################################# Figure 4: Effect of 10% Increase

pdf("Figures/Figure4.pdf",width=8, height=4)
par(mfrow=c(1,2))

data<- read.dta("./Datasets/Other/Swiss_master.dta")
data <- subset(data, ((data$year==2003 & data$cantonid==22) | (data$year==2004 & data$cantonid==10)  | (data$year==2005 & data$cantonid==25)))
data <- subset(data,data$s4_pctforeign > 0)
y <- 100*((data$s4_net_welfare_head) / mean(data$s4_net_welfare_head,na.rm=T))
x <- data$s4_pctforeign*100
general_loess(x,y,1,xlab="Change: % Foreign",ylab="Change in Spending (% of Mean)",c(0,275),"Switzerland",0,.99,trans=1)

data<- read.dta("./Datasets/Other/Swiss_master.dta")
data <- subset(data, ((data$year==2008 & data$cantonid==22) | (data$year==2009 & data$cantonid==10)  | (data$year==2010 & data$cantonid==25)))
data <- subset(data,data$s4_pctforeign > 0)
data <- subset(data,data$s4_pctforeign < oldmax)
y <- 100*((data$s4_net_welfare_head) / mean(data$s4_net_welfare_head,na.rm=T))
x <- data$s4_pctforeign*100
general_loess(x,y,1,xlab="Change: % Foreign",ylab="Change in Spending (% of Mean)",c(0,275),"Switzerland",0,1,trans=0,addplot=oldmax)
legend("topleft", c("Before Reform","After Reform"), lwd=c(1,1),col=c("black","black"),lty=c(3,1),cex=.8)

belgiumdata<- read.dta("./Datasets/Other/Belgium_master.dta")
belgiumdata <- subset(belgiumdata,belgiumdata$year==2003)
data <- belgiumdata
data <- subset(data,data$s5_pctforeign > 0 )
y <- 100*((data$s5_cpas_head)/mean(data$s5_cpas_head,na.rm=T))
x <- data$s5_pctforeign*100
general_loess(x,y,1,xlab="Change: % Foreign",ylab="Change in Spending (% of Mean)",c(0,250),"Belgium",0,.99,trans=1)

belgiumdata<- read.dta("./Datasets/Other/Belgium_master.dta")
belgiumdata <- subset(belgiumdata,belgiumdata$year==2009)
data <- belgiumdata
data <- subset(data,data$s5_pctforeign > 0)
y <- 100*((data$s5_cpas_head)/mean(data$s5_cpas_head,na.rm=T))
x <- data$s5_pctforeign*100
general_loess(x,y,1,xlab="Change: % Foreign",ylab="Change in Spending (% of Mean)",c(0,250),"Belgium",0,1,trans=0,addplot=oldmax)
legend("topleft", c("Before Reform","After Reform"), lwd=c(1,1),col=c("black","black"),lty=c(3,1),cex=.8)
dev.off()


################################## ################################## ################################## 
################################## Figure 5: Immigrant Voter Registration in Belgium

pdf("Figures/Figure5.pdf",width=8, height=4)
par(mfrow=c(1,2))

data <- read.dta("./Datasets/Other/Belgium_master.dta")
data <- subset(data,data$pctnoneunoncit > .01)
data <- subset(data,data$year == 2005)
data <- na.omit(data.frame(data$s2_noneu_enrolled ,data$registered))
x <- (data$data.s2_noneu_enrolled )*100
y <- data$data.registered*100
general_loess2(x,y,.9,xlab="% Change Imm. Beneficiaries, 2003-2005",ylab="Foreign Registration Rate (%)",c(0,45),"",0.025,.975)

### Benefit Usage
data <- read.dta("./Datasets/Other/Belgium_master.dta")
data <- subset(data,data$pctnoneunoncit > .01 & data$benefit_pct > 0)
data <- subset(data,data$year == 2005)
data <- na.omit(data.frame(data$benefit_pct,data$registered))
x <- data$data.benefit_pct*100
y <- data$data.registered*100
general_loess2(x,y,.9,xlab="Proportion Immigrants Receiving Assistance",ylab="Foreign Registration Rate (%)",c(0,35),"",0.025,.975)
dev.off()


################################## ################################## ################################## 
################################## Figure 6: Generosity of cash transfers in Switzerland

pdf("Figures/Figure6.pdf",width=6, height=4)
data <- read.dta("./Datasets/Other/fribourg_aggregate.dta")
plot(data$year,data$ratio,type="l", ylim=c(.6,1),cex.lab=.8,ylab="Ratio",main="",xlab="Immigrant vs Native Per Capita Benefit, by Year",cex.axis=.85)   
points(data$year,data$ratio,cex=.8,pch=18)
abline(v=2004,col="lightgray",lwd=2)
grid()
dev.off()
